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ABSTRACT 


In  this  report,  a  finite  element  method  is  developed  to  approxi¬ 
mately  analyze  umteady  compressible  gas  flow  in  one  spatial  dimension. 
Extension  of  this  analysis  to  three  dimensions  is  possible.  Experimen¬ 
tal  work  to  assist  in  developing  and  verifying  the  model  has  bep>:,i,  but 
will  be  completed  and  reported  at  a  later  date.  The  present  method,  as 
programmed  for  an  jBM  360/65  computer,  is  suited  for  a  variety  of  one¬ 
dimensional  flow  situations,  but  not  for  those  characterized  by  ex¬ 
tremely  large  ratos  of  pressure  change  typical  of  a  weapon.  The  form 
of  the  equations  is  such  that  heat  and  mass  sources,  and  sinks,  fric¬ 
tional  losses,  cross-sectional  area  changes,  and  other  quantities  of 
engineering  significance  can  readily  be  incorporated.  This  analysis  is 
not  intended  to  be  exact,  but  is  to  serve  as  a  mathematical  model  to 
approximate  some  of  the  dynamic  phenomena  in  gas  flow.  This  work  is 
part  of  a  continuing  effort  by  the  Research  Directorate,  Weapons  Labora¬ 
tory  at  Rock  Island,  to  explore  the  basic  mechanisms  of  gas  flow. 
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OBJECTIVE 


The  objective  of  this  study  was  to  develop  a  method  to  approxi¬ 
mately  analyze  unsteady  compressible  flow  by  exploitation  of  certain 
aspects  of  the  philosophy  of  finite-element  techniques  as  applied  to 
mechanical  structures. 

INTRODUCTION 


The  solution  of  the  basic  governing  partial  differential  equations 
of  fluid  mechanics  is  presently  impractical  for  many  important  flow 
situations.  Thus,  approximations  and  idealizations  are  required.  In 
structural  analysis,  a  successful  approach  to  complex  problems  is  to 
first  divide  a  complicated  body  into  elements  that  can  be  individually 
analyzed  and  then  to  account  for  their  interactions.  The  underlying 
purpose  of  the  present  study  is  to  explore  the  extent  to  which  this 
basic  idea  can  be  applied  to  some  fluid  problems. 

In  general,  the  finite-element  concept  is  based  on  representing  a 
continuum  by  a  collection  of  a  finite  number  of  component  parts  joined 
at  prescribed  nodal  points.  Usually  each  component  is  chosen  to  have 
a  relatively  simple  geometric  shape.  Various  field  quantities  are 
locally  described  over  each  element  and  are  assumed  to  be  uniquely  de¬ 
fined  by  their  values  at  the  nodes  of  the  element.  Thus,  each  element 
can  be  considered  to  be  disjoint  for  purposes  of  describing  its  local 
behavior.  Once  the  behavior  of  a  typical  element  is  defined,  a  discrete 
model  of  a  continuum  of  almost  any  shape  with  arbitrary  boundary  and 
initial  conditions  can  be  obtained  by  the  connection  of  elements  in  an 
appropriate  manner.  A  variational  principle  governing  the  problem  at 
hand  is  customarily  formulated.  Such  a  formulation  provides  a  consis¬ 
tent  way  to  generate  analogues  of  the  field  equations,  which  hold  in  an 
average  sense  across  each  element.  However,  the  use  of  a  variational 
approach  is  unnecessary  if  a  well-defined  functional  is  unavailable  for 
the  particular  problem  under  consideration.  Appropriate  forms  of  the 
conservation  laws  can  also  be  usea.  For  viscous  flow,  no  general  varia¬ 
tional  functional  exists,  and  no  attempt  is  made  in  the  present  study  to 
use  a  variational  principle. 

This  study  is  not  designed  to  provide  a  rigorous,  exact  analysis, 
but  rather  an  approximate  model  that  might  provide  some  of  the  charac¬ 
teristics  of  dynamic  gus  flow  at  least  qualitatively.  Some  accuracy  is 
willingly  compromised  for  the  possible  benefits  of  computational  speed 
and  conceptual  simplicity. 

In  many  flow  situations,  series  of  naturally  occurring  element 
types  exist  such  as  bends,  areas  of  surface  roughness,  heat  and  mass 
sources  and  sinks,  and  changes  in  cross-sectional  area.  Seemingly,  the 
equations  governing  individual  elements  could  be  reasonably  joined  in 
some  fashion  analogous  to  that  in  structural  analysis. 

For  the  simple  flow  situation  chosen  to  initiate  this  study,  i.e.. 
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flow  into  a  circular  tube  with  a  piston  blocking  one  end,  the  method  of 
characteristics  would  probably  in  many  respects  be  superior.  Also, 
other  finite-difference  techniques  could  be  applied  to  the  governing 
partial  differential  equations.  However,  many  practical  problems  are 
characterized  by  complex  geometrical  configurations  and  boundary  condi¬ 
tions.  To  reduce  computation  time  by  traditional  methods  requires  non- 
uniform  meshes.  The  setting  up  of  the  difference  equations  can  be 
difficult  in  these  instances.  In  looking  beyond  these  traditional  meth¬ 
ods  and  toward  the  finite-element  approach,  one  can  anticipate  that  a 
versatile  method  could  possibly  emerge  by  which  nonuniform  meshes,  com¬ 
plicated  shapes,  and  complex  boundary  conditions  can  be  handled  in  a 
straightforward  manner. 

APPROACH 


As  a  first  step  in  this  study,  literature  surveys  were  made  that 
included  finite-element  techniques  for  solids  and  fluids,  and  solution 
methods  for  simultaneous  nonlinear  ordinary  differential  equations. 

Very  little  information  was  found  concerning  finite-element  techniques 
applied  to  fluids.  All  the  information  that  was  collected  served  as 
background  and  guidance  in  developing  an  approach.  The  availability  of 
an  experimental  apparatus  to  help  guide  the  effort  from  a  physical 
standpoint  was  considered  to  be  a  valuable  asset.  With  the  use  of  this 
apparatus,  one  should  be  able  to  explore  wide  ranges  in  pressures,  tem¬ 
perature,  densities,  and  velocities.  Although  the  initial  plan  was  to 
study  low  flow  rates,  the  expectation  was  that  eventually  the  analysis 
technique  developed  could  be  extended  to  describe  the  extreme  situations 
encountered  in  weapon  gas  flow.  Thus,  a  test  apparatus  was  constructed 
that  consisted  of  a  highly  instrumented  Ml 6  R’fle  barrel  and  gas  tube. 

Low  flow  rates  in  the  gas  tube  could  be  achieved  by  the  partial  backing 
of  the  gas  port.  Also,  bottled  gas  could  be  used  instead  of  propellant 
gases  from  a  cartridge.  Provisions  were  made  for  simultaneous  measure¬ 
ments  of  temperature  and  pressure  at  many  points  along  the  barrel  and 
gas  tube.  Temperatures  at  various  depths  in  the  barrel  wall  at  several 
locations  could  also  be  monitored.  All  the  test  equipment  has  been  pro¬ 
cured,  and  tests  will  begin  in  the  near  future. 

The  mathematical  approach  in  developing  the  finite-element  method  was 
based  on  the  application  of  macroscopic  balances  of  mass-,  momentum,  and 
energy  to  each  element  and  on  the  requirement  of  continuity  of  various 
thermodynamic  quantities  and  their  derivatives  at  each  node.  No  similar 
approach  could  be  found  in  the  literature. 

In  the  macroscopic  balances,  rate  of  change  of  total  mass,  momentum, 
or  energy  contained  within  a  control  volume  is  equated  to  influx  minus 
efflux  of  that  quantity.  In  the  approach  developed,  each  element  is 
considered  a  separate  control  volume.  Transformation  of  the  macroscopic 
balances  into  such  a  form  that  these  equations  are  expressed  in  terms  of 
pressure,  temperature,  density,  and  velocity  at  the  nodes  (boundaries  of 
the  control  volumes)  is  desirable.  The  right  sides  of  these  equations, 
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which  express  influx  minus  efflux,  are  already  in  this  form;  the  left 
sides  expressing  rates  of  change  of  total  quantities  are  not  yet  in  this 
form,  A  key  problem  is  to  express,  as  accurately  as  possible,  the  mass, 
momentum,  and  energy  contained  within  an  element,  wholly  in  terms  of 
pressure,  temperature,  density,  and  velocity  values  at  the  nodes.  One 
approach  is  first  to  assume  general  forms  for  the  spatial  distributions 
along  each  element  for  any  three  of  the  four  quantities  pressure,  tem¬ 
perature,  density,  and  velocity.  The  fourth  quantity  is  determined  by 
the  equation  of  state.  The  assumed  forms  could,  in  principle,  be  dif¬ 
ferent  for  each  quantity  as  well  as  for  each  element.  One  possible  form 

K  (Pressure  "S 

is  a  power  series  Q.  =  i  C.  (t)xn  where  Q.  =<Temperatureiand  X  is 

\  n=o  1  [Velocity  J 

position  in  one  coordinate  direction.  Other  complete  sets  of  base  vec¬ 
tors  could  also  be  used.  Appropriate  combinations  of  these  forms  are 
then  integrated  over  the  volume  of  an  element  so  as  to  give  total  mass, 
momentum,  and  energy  contained  within  the  element.  The  resulting  ex¬ 
pressions  can  then  be  differentiated  with  respect  to  time,  and  the  left 
sides  of  the  macroscopic  balances  are  then  in  terms  of  pressure,  temper¬ 
ature,  density,  and  velocity  values  at  the  nodes.  Thus,  all  the  balance 
equations  are  now  in  the  desired  form. 


So  that  the  number  of  equations  is  equal  to  the  number  of  unknowns, 
additional  relations  must  be  found.  To  determine  the  values  of  the  un¬ 
known  C.  (t)  requires  the  specification  of  auxiliary  conditions.  Rea¬ 
sonable  Conditions  are  continuity  of  pressure,  temperature,  density,  and 
velocity  at  each  node,  and  continuity  of  first  and  higher  order  position 
derivatives  at  each  node.  A  linear  approximation,  corresponding  to  non¬ 
zero  values  of  C^q  and  only,  would  require  continuity  of  pressure, 


cn 


temperature,  density,  and 
A  parabolic  distribution, 


restrictions  on  derivatives, 
nonzero  values  of 


and  Ci2  only,  requires 


velocity  with  no 
corresponding  to 

in  addition,  continuity  of  first  derivatives. 


Ci.O’ 


Cil 


Continuity  of  higher  order  derivatives  would  be  required  for  higher 
degree  polynomial  approximations  in  order  that  a  sufficient  number  of 
equations  exist.  The  values  of  can  thus  all  be  expressed  in  terms 

of  the  thermodynamic  quantities  at  the  nodes.  These  nodal  values  are 
the  unknowns  in  the  problem. 


For  completion  of  the  formulation,  time  histories  of  the  various 
thermodynamic  quantities  are  required  at  certain  points  in  the  flow 
field.  These  histories  must  be  specified  so  as  to  yield  a  physically 
realistic  problem  and  a  mathematically  determinate  one. 

By  the  approach  outlined  above,  the  governing  partial  differential 
equations  of  fluid  mechanics  are  replaced  by  a  set  of  first  order 
n  o  n  1  i  n  ear,  ordinary  differential  equations  in  a  form  convenient 
for  the  input  of  engineering  data. 
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thrnun  Vh!  Problem  selected  to  develop  the  method  was  simple  flow 
Sr°!SL \f °1  ™  circular  cylinder.  In  addition  to  being  approximately 
h!v,!yEe  °f4.tlow  generated  in  the  test  apparatus,  this  motion  is  also 
characteristic  of  stream  tubes.  The  equations  were  developed  to  account 

LLS13  lLCr0SS'Sec^i0nal  area’  but  initially  a  constant  area  was 
assumed.  Also,  governing  equations  were  developed  where  the  specific 

heat  at  constant  pressure  was  expressed  as  a  third-order  polynomial  in 
temperature,  but  initially  a  constant  value  was  assumed.  The  cylinder 
was  hypothetically  divided  into  a  chain  of  shorter  tubes,  or  elements 
as  shown  in  Figure  1.  The  output  from  one  element  served  as  input  to’an 
IS*?!!!  °ne*  kL  near  distributions  across  each  element  fim  used" 
and  latei  parabolic  distributions  were  assumed.  Pressure  and  tempera- 

leteofetheSfuhpflefl  at  -he  inl6t-  and  ve1ocity  was  specified  at  the  out- 
]®*  tube.  A  °f  sixteen  elements  was  used.  The  qoverninq 

equations  for  this  problem  are  described  in  the  next  section  of  ?h?s  9 


GCM/ERNING  EQUATIONS 

Nomenclature: 

t  =  time  (sec) 

p(x,t)  =  density  (Kg/m3) 

p.j(t)  =  density  at  node  i,  i  =  l,2,--n+l 

v(x,t)  =  velocity  (m/sec) 

v-(t)  *  velocity  at  node  i,  i=l,2,--n+l 

T(x,t)  =  temperature  (°K) 

T - ( t)  =  temperature  at  node  i,  i=l,2,--n+l 
P(x.t)  =  pressure  (newtons/m2) 

P^t)  -  pressure  at  node  i,  i=l,2,— n+1 
L1  =  element  length  (m)  i=l ,2 , — n 

Q(t)  =  rate  heat  energy  enters  an  element  (newton-m/sec)  i=l ,2 — n 
W(t)  =  rate  at  wm'ch  fluid  in  an  element  performs  mechanical  work 
on  its  surroundings 
R  =  gas  constant  (newton-m/Kg-°K) 

a  =  Cp  =  specific  heat  at  constant  pressure  (newton-m/Kg-°K) 

A  =  cross-sectional  area  (m2) 

MT0T  =  total  fluid  mass  in  an  element  (Kg) 

^TOT  “  total  momentum  in  an  element  (Kg-m/sec) 

Ejoj  --  total  energy  in  an  element  (Kg-m2/sec2) 

KT0T  =  total  kl‘netic  energy  in  an  element  (Kg-m2/sec2) 

W  =  mass  flow  rate  (Kg/sec) 

F  =  force  of  the  fluid  on  the  solid  and  is  composed  of  the  sum  of 
all  viscous  and  pressure  forces  (newtons) 
g  =  acceleration  due  to  gravity  (m/sec2) 

<v>  =  average  velocity  of  v  over  a  cross  section 
V  =  v  *■  v:  v  =  time  average  value  of  v  at  a  point 

v3  =  perturbation  of  v  about  7 
=  per  unit  mass 

U  +  PV  =  enthalpy  per  unit  mass 
$  =  potential  energy  per  unit  mass 


A  tube  is  divided  into  n-sections  called  elements  as  shown  in  Figure  1. 
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FIGURE  1 


The  flow  in  each  element  is  assumed  to  be  governed  by  the  following 
macroscopic  balances  for  a  control  volume  (Reference  1). 


"ass:  at  M70T  ■  pi<VAi  -  e1+1<vw>A1+1 


Momentum:  =  (^5T+  Pff)i  -  (^W  +  P*),+1  -  f+  MT0T  g" 


Energy:  £  ET0T  -  [(U  +  PV  +  £  ♦  i)W],  -  C{U  +  PV  ♦  i 


1  <v3: 


+  *)W]i+1  +  Q  -  W 


Where  M 
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pdV 
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L 

UT0T  ‘ 
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p 

0 

0 

C  dT  d  V 
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C  -R 


L  (p2  VR  [(A,  +  3-a2)ti  +  AiTi  +  ^T-j]  +  T2  ~12~  ^1p1  +  ^2P1 


VR 


*  P*] (A.j  +  3 Ag ) D  +  -J2  C~P-jT^ (3Ag  *  A-j)  +  p-|T-j  C 3 A-j  +  Ag ) 3 


Calculations  were  also  made  on  the  assumption  that  Cy  is  a  third-degree 
polynomial  in  T.  (The  results  are  too  lengthy  to  present  here.) 

Now  assume  that  the  cross-sectional  area  is  constant  and  that  the  element 
lengths  are  variable.  The  above-cited  equations,  when  applied  to  the 
whole  system  of  elements,  become,  after  much  algebraic  manipulation,  the 
fol lowing: 


^  p2  =  L-|  {p1Vl  "  P2V2}  “  P1 

[4v£  +  2v-j]p2  +  [4p^  +  2p2]v-j  +  [2p^  +  4p2]v2  = 


(Mass) 

(Momentum) 


12 


-  (4v1  +  2v2}Pl  -  ^  (p2v22  .  p^v^2  +  Rp^  _  Rp^T^ } 


Vt2  v„2  VtV 


r-i  ,  _2__  .  1H2  .  o^R  UT  +  9T  \yn 
^24  8  +  12  12  ^2  2T1^p2 


(Energy) 


Pi  Po  Vp  . 

+  ^vl^T+  ]2^  +  12  ^pl  +  P2^V1 


f  ^12  +  4^  +  12  ^P1  +  P2^V2  + 
Vyr  (2p-|  +  4p2)]t2  = 


8 


Vi 2  v,v 


„  i  — 1 ..  *  JL_  i  1  2  +  o^R  /  „  >  • 

8  24  12  T F  (2T2  41V}pl 

-  {yy  (4P]  +  2p2)}f1 

1  • 

+  r-  Ti  -  p2v2  T2)  +  1  p1v13  -  j  p2v23}  +  ■  ]Al  -1 


For  i  =  2,3 _ n  -  1 


[4v1  +  2Vl^1  +  C4Vl  +  2*i]J1+1 

*  [4p1  +  2p1+13v1  +  [2p.  +  4pf+)^i+l 


'  -  L,  (pi+l  VM  -  ‘V?  ♦  Rp1+1  Tj+,  -  RP.  T,} 

V  •  2  t/2  \i 

[— —  +  —  T  T  *  1  1  Ot-R  /  nj  ,  IT  n1 

L  8  24  12  +  T T  {2Ti+l  +  4VJpi 

+  [ir +  t1  +  --tt1  '  ir  <4Ti.i +  2Tf)]i1+1 


+  rv  (il  +  P-}41)  +  vi+-l  /  +  n* 

Lvi  1  4  12  >  +  T r  (pi  +  pi+l^vi 

*  tvi+l  <T?  *  ^T1)  *  IF  Cp,  +  p1+,)]v, 


*  Clf  ^  +  2p1+i)fli  +  [t#  <2pf  +  4p1+i)]f1+1 


(Mass) 


(Momentum) 


(Energy) 


q  <a(eivi  Ti  -  vi+i  Ti.i> +  ?  vq  -  \  *i+1»  *  ^ 
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Final  ly 


[1]  Pn  +  Cl]  Pn+]  _  j_^  ^nvn  "  pn+i  vn+l^ 


(Mass) 


[4v„  +  2vn+13pn  +  [4vn+,  +  2vn35n+1  +  [4pn  +  2pn+1]vn  (Momentum) 


=  ‘  i2pn  +  4pn+l>Vl  ‘  r  !pn+l  vn+l  '  Vn 


+  Rpn+1  Tn+1  ‘  R  pn  TnJ 


v.2  v2 


rJL.  Vtvf1  +  n  n+^  +  (a~R)  (2T  +  4T  )1p  + 

L  8  24  12  12  u'n+1  4  V *pn 

i/2  w2  VV-.  „ 

Nir  +  ~1T  +  "Tr1  4  7T  (4Tn+l  +  2Tn):1Vl 
+  t”n  (T  4  TT->+  TT  (pn  4  pn+l>^n 
'  CTF  (4pn  4  2pn+l)]tn  4  [TF  (2pn  +  4pn+l)]J+1 

■  -  {v„,i  (w  +  +  w  (pn  4  nn4.i)>v„4.,  + 


n+1  '12 


12  'pn  pn+l;  n+1 


(Energy) 


C  la  (pnvnTn  ‘  pn+lvn+l  Tn+1>  4  \  pnvn3  ‘  I  pn+lvn+l;  4 


On  -  Wn 
^n  n 

AL 


The  previous  equations  are  of  the  form  A(x,t)x  =  F(x,t)  where  A(x,t) 
is  a  3n  x  3n  matrix  function, 

X  =  (vr  P2,  v2>  t2  ••••  pn,  vn,  fn,  Pn+1,  tn+1)T  and  F(x,t)  is  a 

nxl  column  vector  function.  The  initial  state  of  the  gas  is  known  and 
we  specify  it  as,  x  . 


Ai =  4vi +  2vi+i 
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Bi  -  4v1+l  +  2vi 


ci =  4pi +  2pi+l 


=  Zp,  +  4p1+1 


!j_  +  ^1±1  +  Vi  Vit1  +  aii  ( 2T  +  4T  ) 
8  24  12  12  U  H1  V 


V  .  vi+l  .  vi  Vl  .  «-R  +  „  , 

IT  +  —  +  Tc  +  IF  (4Ti+l  +  2Ti> 


„  ,pi  .  pi+l,  .  vi+l  ,  *  , 

vi  (  4  +  12  )  +  12  (pi  pi+r 


Pi  pi+i  vi 

vi+1  %  +  T~^  +  T2  ^pi  +  piH^ 


Ji  =  if  (4pi  +  2pi+l> 


Ki  =  TF  (2pi  +  4pi+l> 


Xi  =  L.  {pivi  '  pi+1  vi+l} 


V“E7<pi+i  vf+l  -pivi2  +  Rpi+l  Ti+1  -RpiTi> 


Zi  "  UT  {a  (piviTi  -  pi+lvi+lTi+l ) 


,  1  ,  1  ,  , 

"2  pivi  "  2  pi+lVi+1}  +  1LT 


n 


Note  that  ,  T-j ,  and  vn+-j  are  specified  functions  of  time.  P-j 

P,  .  h  Ml 

T-j  are  known  and  is  determined  from  =  ^y-  and  p-j  =  ^y--  - 


and 


-  p 


1  1 


1  'P- 


A  total  of  3n  equations  must  be  solved.  The  unknown  parameters  are 

vl»  p2’  v25  T25  p3’  v3’  T3s  ”,pn’  V  TnJ  pn+l  *  Tn+1 '  The  fuictions 
pr  Tl’  vn+l  anc*  their  derivatives  are  given. 

TECHNIQUES  EXPLORED  FOR  SOLVING  THE  EQUATIONS 

A  number  of  approaches  were  taken  to  find  a  stable  algorithm  for 
solving  the  system  of  nonlinear,  first-order,  ordinary  differential 
equations  resulting  from  the  application  of  the  finite-element  method  to 
the  hoTlow  cylinder.  These  are  briefly  reviewed  below: 

1.  The  initial  approach  was  to  replace  derivatives  by  finite  dif¬ 
ferences  and  then  solve  the  resulting  system  of  algebraic  equations  to 
find  pressure,  temperature,  density,  and  velocity  changes  at  each  node 
for  a  small  change  in  time.  Even  with  very  small  time-increments, 
however,  the  solution  failed  to  converge.  A  practical  limit  exists  as 
to  how  small  a  time  increment  can  be.  As  the  increment  is  reduced,  many 
more  calculations  must  be  performed.  Round-off  error  begins  to  grow, 
and  accuracy  is  destroyed  just  as  surely  as  if  very  large  time- 
increments  were  used. 

2.  Next,  a  fourth-order  Runge-Kutta  method  was  tried.  This  tech¬ 
nique  was  stable  for  small  time-rates  of  change  and  for  long  time-spans 
after  flow  initiation.  This  technique  ultimately  became  the  most  suc¬ 
cessful  of  the  numerical  techniques  explored. 

3.  A  modified  Hamming  predictor-corrector  method,  which  is  a  stan¬ 
dard  IBM  integration  routine,  was  also  applied  to  the  equations.  With 
this  technique,  the  integration  interval  was  automatically  subdivided  to 
meet  a  specified  accuracy.  However,  stability  of  this  method  was  not  as 
good  as  that  of  the  fourth-order  Runge-Kutta  algorithm.  Also,  this 
method  was  more  expert ive  to  run. 

4.  A  fifth-order  Runge-Kutta  method  was  also  used.  However,  it 
was  less  stable  than  the  fourth-order  Runge-Kutta  technique.  This 
instability  could  be  attributed  to  the  larger  number  of  derivative  eval¬ 
uations  that  were  required. 
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5,  Various  attempts  were  made  to  linearize  the  system  of  equations 
and  then  to  analytically  solve  a  sequence  of  linear  differential  equa¬ 
tions  with  constant  coefficients  valio  over  one  time-interval.  However, 
this  approach  always  resulted  in  an  eigenvalue  problem;  and  the  corre¬ 
sponding  matrix  is  not  symmetric.  This  approach  apparently  would  have 
been  less  successful  than  the  Runge-Kutta. 

6.  Normally,  pressure  and  temperature  were  specified  functions  of 
time  at  one  end  of  the  tube,  and  velocity  was  specified  at  the  other. 
Various  combinations  of  specified  quantities  at  different  locations  were 
evaluated  for  their  effect  on  stability.  Some  of  these  combinations 
resulted  in  a  coefficient  matrix  that  was  singular.  Two  combinations 
without  this  problem  are  given  below: 

a.  Specifying  the  pressure  and  temperature  at  one  end  of  the 
tube,  and  velocity  at  the  other; 

b.  Specifying  pressure,  temperature,  and  velocity  all  at  one 

end. 


7.  Methods  were  developed  and  applied  to  allow  shock-wave  type  of 
discontinuities  in  the  thermodynamic  quantifies  where  instability  was 
imminent.  This  approach  only  slightly  delayed  numerical  instability. 

8.  Linear  distributions  across  individual  elements  were  replaced 
by  second-degree  polynomial  (parabolic)  distributions  for  purposes  of 
calculating  total  mass,  momentum,  and  energy  within  an  element.  For 
evaluation  of  the  additional  coefficients  C^>  all  first  derivatives 

were  assumed  to  be  continuous  at  each  node.  Computation  time  was  con¬ 
siderably  increased  since  the  coefficient  matrix  was  not  banded.  Sta¬ 
bility  was  not  significantly  increased,  but  oscillations  of  values  at 
the  nodes  were  noticeably  reduced.  In  Appenaix  D,  equations  for  the 
parabolic  distribution  ore  presented. 

9.  A  movable  piston  was  placed  at  the  closed  end  of  the  tube,  and 
the  fluid  velocity  was  constrained  to  be  equal  to  the  piston  velocity. 

The  pressure  acting  on  the  piston  face  was  used  in  Newton's  second  law 
to  determine  this  velocity.  Stability  was  only  slightly  delayel  by  this 
modification.  (Ordinarily,  the  velocity  at  one  end  was  constrained  to 
be  zero.) 

10.  A  C0MC0R  Cl  5000  analog  computer  was  used  to  solve  the  equa¬ 
tions  governing  first  one  element  and  then  two  elements.  Good  agreement 
was  achieved  between  analog  and  digital  results.  With  the  use  of  the 
analog  computer,  the  equations  were  solved  for  quite  high  rates  of  change 
that  led  to  instabilities  in  the  digital  approach.  The  analog  computer 
had  an  insufficient  number  of  multipliers  to  solve  the  equations  for 
more  than  two  elements. 
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11 «.  The  equations  were  made  nondimensional  to  determine  whether 
higher  time  rates  of  change  could  somehow  be  more  easily  accommodated. 
No  appreciable  improvement  was  found. 


RESULTS  AND  CONCLUSIONS 


For  the  assumption  of  linear  spatial  distributions  and  for  the  case 
of  flow  into  an  initially  quiescent  tube  blocked  at  one  end,  the  propaga¬ 
tion  of  a  velocity  wave  and  its  reflection  from  the  closed  end  can  be 
seen  in  Figure  2.  Relatively  small  oscillations  of  pressure,  velocity, 
and  temperature  above  and  below  the  "zero  line"  and  downstream  of  the 
advancing  wave  are  evident  in  Figures  3,  4,  and  5.  These  oscillations 
greatly  diminish  at  nodal  points  when  parabolic  distributions  are 
assumed. 

Of  all  the  numerical  techniques  tried,  the  fourth-order  Runge-Kutta 
algorithm  appeared  to  be  the  most  stable.  However,  even  with  this 
method,  the  solution  became  unstable  if  very  high  rates  cf  change  in  the 
thermodynamic  quantities  were  present.  For  example,  if  the  initial  rate 
of  change  of  pressure  was  about  one  atmosphere  per  millisecond,  the 
velocity  wave  was  propagated  to  the  blocked  end  of  the  tube,  was  re¬ 
flected,  and  was  returned  no  more  than  halfway  before  the  solution  be¬ 
came  unstable.  As  initial  rates  of  change  were  reduced,  the  onset  of 
instability  was  delayed. 

Dependent  upon  how  high  the  rates  of  change  of  thermodynamic 
quantities  are  and  how  long  after  flow  initiation  the  solution  is  de¬ 
sired,  some  practical  problems  can  be  solved  with  the  present  finite-r 
element  method  on  an  IBM  360/55  computer.  The  most  promising  approaches 
for  extending  the  region  of  applicability  of  the  method  lie  in  the  use 
of  a  hybrid  computer  and  higher  order  distributions.  Where  this  method 
is  applicable,  it  is  quite  versatile;  is  well  suited  to  handle  area 
changes,  friction,  and  mass  sources  and  sinks;  and  is  conceptually 
straightforward.  This  method,  moreover,  could  be  extended  to  three 
dimensions. 
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DISTANCE  CM?  txio 


Temperature  versus  position  along  the  tube. 
Lines  beginning  at  high  temperatures 


DISTANCE  W) 


RECOMMENDATIONS  FOR  FUTURE  WORK 


To  further  explore  this  finite-element  method  and  to  extend  its 
range  of  usefulness,  the  following  should  be  undertaken: 

1.  Use  a  large  analog  or  hybrid  computer  to  solve  the  equations. 
Integrations  would  be  more  accurate  than  with  a  digital  computer,  and 
high  rates  of  change  may  be  possible. 

2.  Seek  or  develop  other  numerical  techniques  for  solving  the 
simultaneous  equations. 

3.  Use  a  digital  computer  with  accuracy  greater  than  sixteen 
significant  figures  so  that  small  time-steps  can  be  taken  with  minimum 
round-off  error. 

4.  Apply  the  method  of  characteristics  and  compare  results  with 
the  finite-element  technique. 

5.  Run  experiments  on  the  test  apparatus  developed  and  compare 
results  with  predictions. 

6.  Develop  better  methods  to  approximate  the  total  mass,  momentum, 
and  energy  contained  in  an  element. 

7.  Develop  a  variational  formulation  of  the  problem. 

8  Use  higher  degree  polynomial  distributions  across  each  element. 

9.  Explore  the  use  of  distributions  other  than  polynomial. 
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APPENDIX  A 


Computer  Program  with  Assumed  Linear  Distributions 
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G  LEVcL  d>. 


MAIN 


DATE  =  72IOI 


11/17/Of 


IMPLICIT  RFAL*8!A-H,0-Z) 

C  DEFINE  FILE  9(5C,7A,U,N) 

EXTERNAL  FCT,OGTP 
DOUBLE  PRECISION  L,NUl 

DIMENSION  PI (9) ,RHOI  ( 9} ,NUI ( 9 ) »  T I (9 ) ,X { 2 A ) , C6RY { 24 ) 

DIMENSION  TIMKA)  ,  PRES  I A  )  ,  TEMPE  (  A  )  ,  T  I  M2  I  A  )  ,  PRMT  {  3  ) 

COMMON /BLOCK 1/PI»RHGI ,  T  I »NUI 
COMMON /BLOX/TIMltTI M2, PRES, TEMPE 
COMMON/ B LOC /L, ALPHA, R,QDN,WDN, A 
COMMON /A  SI T/NEPI.NE ,NOUN,NCUPl ,NEM1 
COMMON/OLDE/IUT 

COMMON/DOTE /XDUT1 » XDCT2  , ILNT1 

COMMON/P  T  S/ I P  T  S 

COMMON/SSS/VDOTI.DI SPL 

COMMON /GAS/HG 

COMMON/CCC /CGAM 

COMMON/BBB/DERY 

COMMON/ IN  I TT /ICAMP 

D I SPL=C. CDC 

VDOT1=C.ODO 

IPTS=C 

ILNT  1= 1 

(CAMP=  3 

XDOT 1=C . COO 

XOO  T2  =  C . CDF 

;ut=: 

NkMEL  ist/inpti/ne,a,alpha,r,oelt,niter,wdn,qdk,l 

READ  (b,  INPTl) 

PRINT  1 1  ,NE  * A , ALPHA  »R»DELT»NI TER,  HON, CDN,L 
11  FORMAT! 5X » 1 NE  = ' , I  3 , 5X , » A= ' ,D1 5 . 8 ,5X , • ALPHA= • , 

1F1C.5,/,5X, 'R=' ,F10.5,5X, •DELT=',D15.8,5X, 

2  'NITER=» ,  15, /,5X,fWDN='  ,F10.5,5X,*CDN=  *,F  10.5,  t>X,»L=' ,  F  10.8,  /) 
READ  1 , P I  0 ,RHOI C , XUI 0 , T 1 0 
READ  SOl.NTIM, ( TIM1 (I) , PRES (I) ,l=l,NTIM) 

1  FORMA  I  (  8F 1 C . b ) 

READ  502,MTIM, (TIM2 ( I) , TEMPE ( I ) , I =1 ,  NT  I M ) 

5C1  FORMAT! I b, / , ( 8E 1C. 5 )  ) 

5C2  FORMAT! I  5, /, ! 8F 10. 5 ) ) 

C  HG= 2. C*3. K*o SORT! 7.9172  90-6) *A 18 .6/ ( .002*DSGRT ( 3. 1  ADO  )  ) 

C  HG=-HG 

HG=  0 . CDC 

CGAM=ALPHA/( ALPHA  -  R) 

PRINT  779, HG 

779  FORMAT!  '  HG=',D12.5) 

NEP 1=NE  +  1 
DO  5CC  I  =  1 , Nt P 1 
PI! I )=P!C 
NUI  (  I )  =  XUlr 
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G  LEVFL  2,;. 


MA  f  N 


5CC 


C  MtW 

C 

C 


5C 

C 

c 

c 

c 

c  50 

C  NEVv 
51 

C 

C 

C 

C 

C 

C 

C 

C 


522 


r  n i j = rr  c 

I(  F  ) "  P  I  C'  /  (  R  *  7 1  r  ) 

M0UN=3*NE 
N0Uei=N0lA'  +  1 

N  E  M 1  =  N  E  _  t 
BOUNDARY  C0NDI7IDNS 

temporary 

READ  1 , p  j 
BEAD  1,  RHOI 
READ  1,N0I 
READ  1 f  T  j 

temporary 

X(  U=NUI(  1) 

X ( NOUN  -  l )  =RHOi FNEP1 ) 

X(  NOUN)  =  n  (NEP1 ) 

IF(NE.EO.l)  GO  70  51 
DO  50  F=2f NE 
ITE=3*i 

X(ITE  -  4 ) =RHO I ( I ) 

2!JI£  '  3)  -NUKI) 

X(  I  TE  -  2  )  =  7F  (  I  ) 

DO  50  I  =  2 1 NE P 1 
F  TE=  3*  F 

XUTh  -  5  )=rho  I(J) 

*  ‘  1 F  -  4 )=NU: ( J ) 

X<ITE  -  3 ) =71 | I ) 

;?^OARY  CONDI  7 IONS 
READ  1 , PRM7 

CALL  RUNGEIPOMT.X.OeRY.KCU.FCT.OI 

WRITE  (O' I) PI 

1  =  2 

WRITE  ( 9*  F ) RHOI 
1  =  3 

WRI7EF  9' F ) NUI 
1  =  4 

WRITE (9* I  )tj 
PRIN7  522.FP7S 
F0RMA  7F.  IP  TS= ' , I 3) 

CALL  Exit 
END 


DATE  =  72101 


11/17/00 
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ooooooo  oooo 


LEVEL  20 


OUTP 


CATE  =  72101 


11/17/00 


SUBROUTINE  OLTP! T,X,DERY,NDI M,PRMT) 

IMPLICIT  REAL*8!A-H,0-Z) 

DOUBLE  PRECISION  NUI,L 

01 MEN  SI ON  X  (  24  ) ,DERY<24) »PRMT(3) 

DIMENSION  PI (9) ,RH0I (9) ,NUI ( 9 ) , TI (9 ) 

COMMON/ASIT/NEPltNE  »NOUN» N0UP1 ,  NEMl 
COMMON/BLOC  /L  , ALPHA , R  ,  QON  ,  WON  ,  A 
COMMON/BLOCK  1/P  I ,RH0I ,71  ,NLI 
CUMMPN/OLOE / 1 UT 
COMMON /P  TS/ 1 P  T  S 
COMMON/ SSS/VDOT 1,01 SPL 
IFOWDUUT,  5).NE.O)  GO  TC  401 
NUI ( 1 ) -  X  (  1 ) 

T  I  (NEP1 )  =  X( NOUN ) 

C  NEW  BUUNDAR Y  CONDITIONS 

RHfll  (  NEP1  )=X( NOUN  -  1) 

C  FOR  INTEGRATING  END  VELOCITY 

STORE=P I ( N£P 1 ) 

P I ( NFP 1 ) =R#RHOI ( NEP 1 ) *  T I (NEP1 ) 

Vt)OTS=VOCTl 

IFIT.GT.l.OD  -3) V00T1  =  VD0T1  +  PRMT 1 3 > *  I  I P I t NEP 1 )  +  STORE)/2.DO 
1  -  101000.00*7.417290-2 
0 1  SPL  =  0  I  SPL  +  PRMTI  3)*!  VDCTS  +  VDCTD/2. 

FOR  INTEGRATING  END  VELOCITY 
IFINE.bO.l)  GO  TC)  61 
00  51*  1  =  2, ME 
ITC=3*I 

R HO  1  (  1  )  =  X(  1  Tt  -  4) 

NUI !  I  )  =  X( ITE  -  3) 

TI (  1  )  =  X( ITE  -  2) 

64  PI ( I)=R*RHOI C I ) *TI ( I  ) 

00  54  I=2»NEP1 
I TE=3* I 

RHOI ( I  )  =  X( I TE  -  5) 

NUI ( I )  =  X(  I  7E  -  4) 

Tit  I )  =  X ( ITE  -  3) 

54  PI ( I )=R*RHOI ( I )*TI (I) 

NEW  BOUNDARY  CONDITIONS 
61  PRINT4), T,L,DISPL,UDN, WON 

9  FORMAT!'  T=  «,  F12.8,  *L=  ’  ,F12.8,  '  DISPL=  D12.5, 

1  '  QDN=  ' ,D 1 2 . 5 , '  WDN=  ',  D12.5,/) 

IPTS= IPT  S  +  1 
C  IF! IP  TS.GT.50)  GO  TC  709 

C  WRITE ( 9' IPTS) T,PI »RHOI »NUI  ,TI 

709  PRINT  3 

3  FORMAT!'  N0DE',5X,'  PRESS.  0IST.',4X,'  DENSITY  DIST.', 

1  5X , '  VELOC.  OIST. '  ,6X,'  TEMP.  D I  ST . *  , / , 1 3X  , 

2  'NT/M**2' ,11X, »KG/M**3',12X, 'M/SEC' , 1 3X , 'KELVIN'  ) 
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FUNCTION 

1  F$< I,J,RH01,NUl,Tl,RH02,NU2,T2) 

IMPLICIT  REALMS ( A-H , C-Z ) 

DUUBLE  PRECISION  NU1,NU2,L 

COMMON/BLOC /L, ALPHA, R,QDN, WON, A 
COMMON/DOTE/XDOT1 ,XD0T2  ,  ILNT1 
GO  TO  (1,2,31,1 

1  GO  TO  ( 11,12,13,14,15,16,171  ,J 

2  GU  TO  (21,22,23,24,25,26,27) ,J 

3  GO  TO  (  31,32,33,34,35,36,37)  ,J 

11  FS=1.DC*L 
RETURN 

12  FS=0.000 
RETURN 

12  FS=O.ODO 
RETURN 

14  FS= 1.0C*L 
RETURN 

15  FS=0. COO 
RtTURN 

16  FS=O.COC 
RETURN 

17  FS=2.4( RH01*NU1  -  RH02*NU2 
1-  XD0T2*RH02  +  X00T1*RH01> 

RETURN 

21  FS=(4.*NU1  +  2.*NU2)*L 
RETURN 

2?  F S= ( 4 . *RHO 1 ♦  2.*RH02)*L 
RETURN 

23  F S=  0 . COO 
RETURN 

24  FS=<4.*NU2+  2.*NU1)*L 
RETURN 

26  FS“ ( 2 .*RH01+  4.*RH02)*L 
RETURN 

26  FS=0. CDO 
Re  TURN 

27  FS=  -12.*(RH02*NU2**2  -  RH01*NU1**2  +  R*RHC2*T2  -  R*RH01*T1) 
1  -( X00T2*NU2*RH02  -  XD0T1 *NU1 *RH01 1 *12. 

RETURN 

31  FS=(NUl*NUl/8.+  NU2*NU2/24.  +  NU1*NU2/12. 

1  + ( ALPHA  -  R)*(2.*T2  +  4.  *T1 1 / 1 2 , 1 *L 

RETURN 

32  FS=(NUl*(RH01/4.+  RH02/12.)  +  NU2*( RH01/ 12. +  RWC2/ 12. ) )*L 
RETURN 

33  FS= 1 ( ALPHA-  R ) * { 4 . *RHO 1  +  2 . 4RH02 1 / 1 2 . 1  *L 
RETURN 

34  FS=(NUl*NUl/24.-f  NU2*NU2/8.  +  NU1*NU2/12. 
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1  +(  ALPHA  -  R)*(4.*T2  +  2.  *T1  ) /12.  )  *L 
RETURN 

35  FS=(NU2*(RH01/12.+  RH02/4.)  +  NU1*  ( RHC 1/ 12 . +  RHC2/12.)!*L 
RETURN 

36  FS~  (  ( ALPHA-  R  ) * ( 2. *RH01  +  4. *RHG2 ) / 12 . ) *L 
RETURN 

37  FS= ( ALPHA*! RH01*NUt *T l  -  RH02*NU2*T2> 

1  +  RH01*NUl**3/2.  -  RH02*NU2**3/2. >  +  (QDN  -  WDN)/A 

2  -  ( X00T2*RH02*(NU2*NU2/2.  +( ALPHA  -  R)*T2) 

3  -  XDQTl*RH01*(NUl*NUl/2.  +  (ALPHA  -  R)*TD) 

RETURN 


¥ 
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G  LEVEL  20  GAUSSE  DATE  =  72101  11/17/00 

SUBROUTINE  GAUSSE  (X,NJ 
IMPLICIT  REAL*8 ( A-H ,0-Z ) 

CGMMON/ASTHT/  A 

C  USES  GAUSS  ELIMINATION  METHOD  TO  FIND  SOLUTION  TO  SYSTEM  OF  EQUATIONS 

DIMENSION  X(24) ,A(24,25) 

NM 1=N  -  1 
NP 1=N  +  1 
DO  1  K=1,NM1 

C  FIND  FIRST  ROW  WITH  A  BIG  ENOUGH  ELEMENT 

KP5=K  ♦  5 
00  ICC  I  I=K»KP5 

I F {  DABS ( A ( 1 1  ,K  )  )  .GT.  1.0D-30)  GO  TO  101 
ICC  CONTINUE 
PRINT  99 

99  FORMAT!  '  ALL  ELEMENTS  ARE  PRACTICALLY  2ER0  IN  MATINV'  ) 

CALL  EXIT 

101  IF  (II.EQ.K)  GO  TO  103 
C  MAKE  ROW  INTERCHANGE 

KP7=KP5 

N£ND=  MINC(KP7,NP1) 

DO  1C4  JJ=K,NEND 
8= A  (  K  ,  J  J  ) 

AIK, JJ)=A! II, JJ) 

A! 1 1 , JJ  J=B 
104  CONTINUE 

IF!  NEND.EQ.NP1)  GO  TO  103 
8= A ( K , NP 1  ) 

A ( K , NP  1  )  = A ( 1 1 ,NP1  ) 

A(JI,NP1  )=B 

C  INTERCHANGE  COMPLETED 

103  B=A { K , K  I 

DO  2  J=K,NEND 
A(K,J)=A(K,J)/B 
2  CONTINUE 

IF  ( NEND.EQ.NP1 )  GO  TO  106 
A ( K , NP1  ) =A ( K  » NP1  )/B 
1C6  KP 1=  K+l 

NST0P=KP1  ♦  3  -  MOD ( KP1 ,3 ) 

I F ! NSTOP  ,GT.  N)  NSTOP=N 
DO  1  I=KP1, NSTOP 
B= A  (  I  ,  K  ) 

DO  17  J=K  , NEND 
17  A! I , J )= A ( I,J)  -  B*A!K,J) 

IF  ( NEND.EQ.NP1 )  GO  TO  1 
A ( I ,NP  1 ) =A ( I ,NP1 )  -  B*A !  K ,NPl ) 

1  CONTINUE 

107  X(N)=A(N?NP1)/A(N,N) 

DO  3  K=1,NM1 
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KP 1=  N  -  K  +  1 
X ( N-K )  =  A { N  -  Kf NP1 ) 

NSTQP=  KP1  +  M0D( K » 3 )  +  3 
M S TO P  =  M I  NO { NSTOP  T N ) 

DO  3  J=KP1, NSTOP 

3  X( N-K )=X( N-K )  -  A ( N  -  KtJ)*X(J) 
KCTUKN 
END 
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SUBROUTINE  UPDA TE I T , X , Y ,F  ,FD ) 
IMPLICIT  REAL*8(A-H,0-Z) 
DIMENSION  X ( 4 ) *  Y ( A ) 

DU  1  1=1,4 

IF( T  .LE.  XI  I)  )  GO  TO  2 

1  CONTINUE 
1  =  4 

2  1=  1  -  1 
IFU.EG.C)  1  =  1 

F=  Ym  +  (T-xm)*(y(i  + 

FD=  (  Y (  I  +  1)  -  Y<m/(XH  +  l)  - 
RETURN 


1)  -  Y< I ) )/ (X ( I 
Mil) 


+ 


1) 


XII)  ) 
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SUBROUTINE  FCTI T,X,XD0T) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  X(24) ,XD0T(24) ,PI (9) ,RHCI (9) ,NUI (9),TI (9) 

0IMENSI0N  A ( 24*25) 

DIMENSION  DER Y ( 24 ) 

DOUSLE  PRECISION  NUI ,L 
C0MM0N/BL0CK1/PI ,RHQI ,TI ,MUI 
1/BL0X/TIM1,TIM2,PRES,TEMPE 
COMMON/BLOC /L, ALPHA, R,QDN, WON, DUM1 
2/AS  IT/ME? 1,NE, NOUN, N0UP1»NEM1 
CL'MMON/OOTE/XDOTl  ,XD0T2  ,  ILNT1 
DIMENSION  PRE  S ( 4 ) , T I  Ml ( 4 ) , T EMPE (4 ) , T I  M2 (4) 

COMMON  /ASTHT/A 
COMMON/ SSS/VDOT 1,0 1 SPL 
COMMON/GAS/HG 
COMMDN/CCC/CGAM 
C0MM0N/8BB/0tRY 
COMMON/ IN  I T  T/ICAMP 
ICAMP=  ICAMP  4-  1 
N=NOUN 

C  DO  8888  1  =  1, N 
C  00  8888  J= 1 , NOUP 1 
C8880  A ( I , J )=O.CDC 
C  NEW  BUUNOARY  CONDITIONS 
VM=O.COC 

NEW  OOUNOARY  CONDITIONS 
FOR  TRAVELING  WAVE 

IF  ILNT1  IS  GREATER  THAN  ZERC  THEN  NC  TRAVELING  WAVE 
I F ( ILNT1.GT.0)  GO  TO  102 
XNE=NE 
XOO  T 1  =  0 . CDC 
X  DO  T2=37o. O'*  8DQ 
L=T*XD()  i'2/XNF  J-,CC5 
IF(L.GE..C8DG)  GU  TO  101 
GO  TO  102 

101  XU0T2=C . COO 
ILNT  1=1 

FOR  A  TRAVELING  WAVE 

102  CALL  UPDATE! T,TIM1, PRES, PI (1) , PD) 

CALL  UPDATE (T,TIM2,TEMPE,TI(1) ,TD) 

THESE  EQUATIONS  ARE  FOR  THF  ASSUMPTION  THAT  THE  FLOW  CN  THF  TOP  ELFMENT 
IS  ISENTROPIC 
rati=cgam 

PI  (  1 )  =  ( T I (i)/288. 16)**(RATI/(RATI  -  l.))*10100C. 

PD=101C00.*RATI*(T{ ( 1 ) /Z88. 16) *#( 1. / ( RAT  I  -  1.)) 

1  *TD/( 288. 16*(RATl  -  1.)) 

THESE  EQUATIONS  ARE  FOR  THE  ASSUMPTION  THAT  THE  FLOW  ON  THE  TDP  ELEMENT 
IS  ISENTROPIC 
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RH0I(  l)=Pl(l)/(R*TI  (D) 

RH00=PD/(R*TII 1) )  -  PI(1)*TD/(R*TI(1)*TI(1)) 

C  NEW  BOUNDARY  CONDITIONS 
C  NUI ( 1 >=266011. 3DC*T 

NUI (  1  )  =  X( 1 ) 

T  I  (NEP1  )  =  X(N) 

COMPUTATION  FOR  SHOCK 
SHOCK  INSERT  A 
COMPUTATION  FOR  SHOCK 
542  VDDT=0.CDC 

0DN=L*HG*( (Till)  +  TI<2>)/2.  -  286.16) 

IF(QON.GT.O.CDO)  Q0N=0.0DC 
GO  TO  776 

IF  (T.LT.l.D-3)  GO  TO  775 
NU I ( NEP 1 ) =VD0T1 
VDO 7=7.917291 70-2* P I  I NE PI ) 

775  RHO I ( NEP  1 )  =  X ( N  -  1) 

IF(NE„CG.l)VM=VDOT 
IF(NE.EO.l)  GO  TO  27 
DU  U  1=2, NE 
I TE=  3* I 

RHO I ( I )=X(I TE  -  4) 

NUI ( I )=X( ITF  -  3) 

11  T I ( I)=X( ITE  -  2) 

DO  11  I  =  2 , NEP  1 
175=3*1 

RHOK  I)  =  Xt  ITE  ~  5) 

NUI ( I !=X( ITE  -  4) 

11  Till  )  =  X ( ITE  -  3) 

SET  UP  FIRST  3  EQUATIONS 
27  DO  1  1=1,3 
JJ=1 

DO  2  J=4 , 6 
C  JJ=J  -  3 

IF  l  NE  .NE.  1 ) J J  =  J  -2 
IFINE.EO.l  .AND.  J.EQ.5)  GO  TC  2 
IFCNE.EO.l)  JJ=JJ+  1 

A(  I, JJ  )=FS(I*J, RHO 1(1), NUI() ) , T I ( 1 ) ,RHOI 1 2 ) , NU I ( 2 ) , T 1 1  2 )  ) 

2  CONTINUE 

C  IF(NE.EQ.l)  GO  TO  771 

NEND=MINC<  8,N) 

DO  3  .)P  =  5,NEND 

3  A I  I , JP )  =  C . CDC 

C  NEW  BOUNDARY  CONDITIONS 

7  71  A I  I tNOUP 1 ) =r  S ! 1,7, RHO I (l)  ,NUI (1)  ,  T I  ( 1 ) , RHO I  1 2  ) , NU  I  ( 2  ) , T  I  ( 2  )  ) 

1  -  RHOD*FS( I ,l,RHOi  ( 1) »NUI ( l> ,TI (1 ) ,RHCI (2), NUI ( 2 ) , T I (2)  ) 

2  -TDYFSI 1 , 3, RHO I (1),NUI(1),TI(1) ,RHOI ( 2 ) , NU I ( 2  > , T I ( 2 ) ) 

3  -  VM*FS(I,2,RH0I(1),NU III), Till), RHO l (2), NU 1(2), 11(2)} 
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A(  i  ,1  )=FS(  1 ,2»RHGI (  1)  ,NUI  (1.)  ,TI  (1 J ,  RHO I C  2  )  ,NUH  2 )  ,T  II 2 )  ) 

1  CONTINUE 

FOR  A  SHOCK  CONDITION 
INSERT  B 

FOR  A  SHOCK  CONDITION 
FOR  A  TRAVELING  WAVE 
XDOT 1  =  X  DU  T  2 
FOR  A  TRAVELING  WAVE 

FIRST  3  EQUATIONS  ARE  SET  UP 
40 CC  IF(Nt.EC.l)  GO  TO  7C0 
IFINE.EQ.2)  GO  TO  600 
DO  50  11  =  2, NEM1 
COMPUTATION  FOR  SHOCK 
INSERT  C 

END  COMPUTATION  FOR  SHOCK 

668  QDN=L*HG*I  ( T I  (  II )  +  THII  +  l:)/2.  -  288.16) 

IF(QDN.GT.O.ODO)  Q0N=0.0D0 
ITE=3*( I  I  -  1) 

I TE 1=  3* (  1 1  -  2) 

I  TE 1=  3*(  I  I  -  2)  +  1 
NEW  BOUNDARY  CONDITIONS 
DU  51  1=1,3 
DO  52  J=  1 ,6 

52  A (  !  TE  +  I ,  I  TE  1  +  J)  =FS(  I  »  J»RHOI  ( 1 1 )  ,Nl)I  (  1 1 )  ,T  I  (  1 1 ) » 

1  RHOIIII  +  1 )  ,NUI  > I  I  +  1 )  ,  T I (  1 1  +  1)) 

NEND=  4 

I F (  1 1  .C0.NFM1)  NFND  =  3 
DU  53  J= 1 , NEND 

53  A (  I  TE  +  I ,  l TF 1  +  6  +  J)=O.ODO 

51  A(  I  TF  +I,N0UPl)=FS(I,7,RHCnm,NUI(II),THm, 

1  RHOIIII  +1 ) ,NUI ( II  +  1 )  ,  TI I  1 1  +1)) 

COMPUTATION  FOR  SHGCK 
INSERT  D 

COMPUTATION  FOR  SHOCK 
5C  CONTINUE 

NEW  RGUNOARY  CONDITIONS 
GO  TO  7CC 

NEW  BOUMOARY  CONDITIONS 

THIS  IS  USED  WHEN  BOLT  OR  PISTON  IS  ALLOWED  TO  MOVE 
X$TORE=L 
XD0T2=NU I ( NEP 1 ) 

L=L  +  DISPL 

WDN= ( P I ( NEP 1 )  -  101CCG.DG)*DUM1*NUI (NEP1) 

THE  ABOVE  IS  USED  WHEN  THE  BGLT  OR  PISTON  ARE  ALLOWED  TO  MOVE. 
SET  UP  THE  LAST  EQUATION 
QDN=  L*HG*(  (  TUNE)  +  T!(NEPl>>/2.  -  288.16) 
lF(QDN.GT.C.ODC)  QDN=O.ODC 
600  DO  60  1=1,3 
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INT  =  o 

DO  61  J=l,6 

C  NEW  BOUNDARY  CONDITIONS 
C  IF(J.Eti.4)  GO  TO  61 
IFIJ.EC.5)  GO  TO  61 
C  NEW  BOUNDARY  CONDITIONS 
INT= INT  +  1 

AIN  -  3  +  I,N  -  5  +  INT)=FS(I »  J  »RHOI (NE) » NU I (NE)»TI (NE)» 
l  RH0HNEP1)  ,NUI  (NFP1)  ,TI INEP1)  > 

61  CONTINUE 

C  NEW  BOUNDARY  CONDITIONS 

AIN-  3+  I»N0UP1)=FS(I,7*RHQI (NE)»NUI (NE)»TI(NE),RHOI(NEPl)t 

1  NU 1 1  NFP  1  >  t  TI INEP1 )  )  -  VDCT*FS(I,5,RHCIINE),NUI(NE),TKNE), 

2  RHOI (NEPl) fNUI (NEP1) »  T I (KEPI ) ) 

C  1  NUI I NEP 1 ) *  TI I NEPl ) ) 

C  AIN  -  3  +  I » N  )  =A I N  -  3  +  I  »N)  *  RH02*FS(I ,4,RHCI (NE) »NUI (NF), 
C  1  TIINE), RHOI (NEF1) , NUI (NEPl) ,TI (NEPl) ) 

C  NEW  BOUNDARY  CONDITIONS 
60  CONTINUE 

C  USED  WHEN  PISTON  IS  ALLOWED  TO  MOVE 
C  XOO  T2=C . COO 

C  L=XSTORE 

C  WDN=0 .0 

C  USED  WHEN  PISTON  IS  ALLOWED  TO  MOVE 
C  IF! MODI ICAMP , 4 ) . NF . C )  GO  TO  7G0 

C  DIMENSION  STT I  24 ,25 ) 

C  DU  50 5 T  1  =  1, N 

C  DO  5050  J= 1 ,NOUP 1 

C5C5C  STT I  I , J ) - A I  I , J ) 

7C0  CALL  GAUSSE I XOOT  »N) 

C  I F ( MODI ICAMP,4).NE.C)  GO  TO  7000 

C  DO  5051  1=1, N 

C  DEV=0 . CDO 

C  DU  5052  K= 1 , N 

CSC 52  Ot V=DEV  +  STT I  I , K )«XDOT I K > 

C  ERRQR=(OEV  -  STTI I ,N0UP1 ) ) *100. /DEV 

C  PRINT  9040, DFV, ERROR, XDOTI I) 

C904C  FORMAT!'  LEF T=  •  ,01 2 . 5 , '  E RROR=  ' , D12 . 5 , '  XDOTII)  =',D12.5) 
C5C51  CONTINUE 
70CC  CONTINUE 
RETURN 
END 
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SUBROUTINE  RUNGE ( PRM T , X , DER Y , N. FC T » CuTP ) 

IMPLICIT  REAl*8U-H,0-Z) 

01  MEN  SI  ON  X  (  2-4  >  ,DERY(24)  ,AK1  (24)  ,AK2!24)  ,  AK3  1 24  } ,  AK4  (  24  )  ,  STORE  (  24  ) 
1  * PRM  T ( 3 )  *  AK5 ( 24 } f AK6 I  24  I 

t=prm  r 1 1 » 

DELT=PRMT<  3) 

SUM  =  0ELT<‘.50C 
FINAL=PRMT( 2) 

ICC C  CALL  FCT! T,X,D£RY) 

CALL  OUTP! T,X,OERY,N,PRMT  ) 

IF( T.GE. FINAL)  RETURN 
00  1  I  =  1  *  N 
AKK  I  )  -=0bR  Y{  I  )  *DELT 
TEMP=  X (  I  ) 

STORE! I )  =  TEMP 
1  X<  I  )=TEMP  <•  AK1!  I  )*.5D0 
T=T  +  SUM 

CALL  FC,T(  T»X»OEKY) 

Dii  3  I*  1,N 

AK 2 ( I  )  =  DELT*DERY( I ) 

3  X( I )  =  A K 2 ( I )*.5DC  ♦  STORE! I ) 

CALL  FCT! T  »  X  »OERY ) 

OU  5  I  =  UN 

AK 3 ( I  )  =  D£LT*DERY( I ) 

6  X! I )  =  S  T  0  R  E  (  I  )  +  AK3! I  ) 

I  =  T  SUM 

CALL  FCT(T,X,OERY) 

00  7  I  =  1,N 

AK  4 ( I )=DELT*DERYt I ) 

7  X ( 1 )  =  S TORE ( I )  +  (AKl(I)  +  2.00*(AK2!l)  ♦  AK3(I))  +  AK4(l)>* 

1  .1666666666666666700 

GO  TO  1CCC 
END 
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Existence  and  Uniqueness  Proof 


Because  of  the  difficulty  in  finding  stable  solutions  of  the  simul¬ 
taneous  equations  in  this  study,  some  work  was  done  to  establish  existence 
and  uniqueness  of  solution.  The  arguments  presented  below  demonstrate 
that  for  a  single  element  and  for  very  specific  end  conditions,  a  closed 
interval  exists  about  the  initial  time  for  which  a  unique  solution  is 
present.  No  attempt  at  generalizing  the  argument  to  n-elements  and  very 
general  end  conditions  is  made. 

First,  a  general  theorem  from  Reference  3  is  used. 

Defn:  |x|  =  ^  |xi| 

where  x  is  an  n-dimensional  vector. 

Defn:  f(t,x)  satisfies  a  Lipschitz  condition  on  a  domain  D  of  a 
(t,x)  space  if  and  only  if  a  K  >  0  exists  such  that  | f (t ,x-j )  -  fU^)! 

1  K  |x^  -  x^j  for  each  (t,x^)  and  (t,x?)  in  D. 

Gi ven : 

(E)  x  =  f(t,x) 
x(t)  =  C 


Thrm:  Suppose  f  e  (C,  Lip) 

(i.e.,  continuous  in  t  and  Lipschitz  in  x) 
on  the  rectangle, 

R:  |t  -  t|  -  a,  |x  -  c|  -  b  (a.b  >  0) 
and  let  M  =  max | f ( t ,x) |  (t,x)  e  R 

then  a  unique  solution  exists  ip  e  C1  (C1  is  the  set  of  all 
functions  haveing  one  continuous  derivative)  of  (E)  on  |t  -  t|  _  a  such 
that 

'Pi  t)  =  f(t,i|») 

*( t)  =  K 


and 

a 


min(a, 


i) 
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For  the  domain,  D,  of  the  existence,  one  can  show  that  if  the 
partials  with  respect  to  x  are  continuous  for  each  (t,x)  e  D,  then  a 
Lipschitz  condition  holds  in  D. 

One  then  applies  this  condition  rather  than  show  a  Lipschitz  con¬ 
dition  directly. 


Let 


By  the  mean  value  theorem,  if  —  exists  componentwise  for  each 

dX 

x*  e  B  -  {x1  :  x'  =  (1  -  x)x  +  Ax,  0  <  a  S  1}, 
an  x*  e  B  exists  such  that 

f(t,x)  -  f(t,x)  =  (f|)x=x*  (x  -  x) 


NOTE:  The  domain  of  consideration  in  the  existence  theorem  is  convex; 
therefore,  one  can  apply  the  mean  value  theorem. 


Let 


M*  =  max 


max 


i  >J  (t,x)eD 


,!!i, 

3X  • 


3f-  3f, 

Since  — —is  continuous,  Itch  is  implied  to  be  continuous;  and  D  is 
3xi  <^x. 

closed  and  bounded.  Therefore,  M*  exists. 

| f ( t , x )  -  f(t,x)|  =  S  | f ,•  ( t,x)  -  fi( t,x) I 

J=1  0  J 


n  n  3f  •  ' 

=  £  1 2  (^-)  (xi  -  x . ) ! 

J=1  i=l  8xi  x=x*  1  1 
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"  "  K^r)  l 

j=l  i=l  dXi  x=x* 


lxi 


n 

<nH*E 

i=l 


n  M*  jx  -  x 


A  K,  namely  K  =  n  M*,  exists. 

Now,  consider  the  gas  equations,  and  put  them  in  a  form  so  that  the 
previous  results  can  be  applied. 

The  end  conditions  are 

p-j  ~  a-j  b^jt  ,  T-j  -  c^  »  v2  *•  0 

a^  ,  b-j  ,  Cj  are  constants. 

After  substitution,  the  differential  equations  to  be  solved  are 

(1)  p2  =  -b]  +  £  (a1  +  b1t)v1 

(2)  v]  =  -  {2b]v1  +  (a1  +  b]t)v12 

+  T  Hai  +b1t)v12  +  Rp272  -  R(a1  +  b^jT^}  /  (4P]  +  2p2) 

2 

(3)  [74-  +  ^4T2  +  2T1  +  +  lf^vlvl 

2 

+  (2p,  +  4o2)]t2  ‘  -  t-g-  +  (2T2  +  4T,)]b, 

4  r  +  \  pVi3^ 

By  forward  substitution,  equations  of  the  following  form  result: 
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f>2  ~  (i*»  P25  »  Tg) 

v-j  =  p£»  ^2^ 

^2  ~  ^3  p2’  ^1*  "^2^ 


Initial  conditions  are 


p2  '  p20 
v,  (0)  =  v1Q 

t2(0)  =  t20 


The  only  points  where  the  parti als  do  not  exist  and  are  not  contin¬ 
uous  are  at  p2  =  -2p^  and  p^  =  -2p2.  Since  F-j ,  F2,  and  F3  are  continu¬ 
ous  for  all  t,  one  can  choose  a  >  0  to  be  any  value.  Also, 


p-i(O) 

P2(0j  f  — 2 —  or  p2(0)  f  -  2p -j (0) 


(The  conditions  above  would  be  physically  unrealistic) 

which  implies  that  a  b  >  0  exists  such  that  |p2  -  p2q|-  +  |v-j  -  v^j 


+  |T2  -  T20|  5  b  and  the  parti als  of  F^,  F2>  and  F3  are  finite. 


Existence  and  uniqueness  are  present,  at  least  on  the  interval, 
[0,  a],  where 

•  /  b\ 
a  =  rmn(a,pj) 

and  M  =  max(  jF^  +  |F2|  +  |F3|) 
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M  is  determined  on  the  domain: 


\t\  <  a 

I p2  ‘  p2ol  +  IV1  "  v1ol  +  lT2  "  T2qI  '  b 
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APPENDIX  C 

Analog  Computer  Block  Diagram  for  a  Single  Element 


44 


_  "  .<  -  -  .  ~  ?’ 

C0R  0ME  rlfMCMT 

DIFrEftCHtlKV  M>Aiy*£R  SV.OC.H. 

m 


In  the  following  discussion  superscripts  denote  different  functions 
and  not  derivatives.  Define  p  for  the  first  element  by  the  following: 


(x)  =  p1  +  L-P-1-  -  c1  L)x  +  c^2 


Note  that: 


p(x)  (0)  =  p.j  and  p ( 1 ) (L)  =  Pg 

Rather  than  letting  c^  =  0  and  having  a  straight-line  approximation 
for  the  density  on  the  first  element,  one  may  choose  c-j  such  that,  if 
one  extended  the  definition  of  p  over  the  second  element,  then  p  would 
have  the  value  p3  at  x  =  2L  (i.e.,  p( 1 ) (2L)  =  p3). 

The  appropriate  value  for  c^  is  given  by: 
c-|  =  (p3  -  2f>2  +  p-j)  /  2L2 


Choosing  c^  as  shown  above  yields  better  results  for  the  derivative 

between  elements  (1)  and  (2)  since  the  curve  is  now  bent  in  the  proper 
direction.  For  the  remaining  elements,  the  derivatives  will  be  forced 
to  be  continuous  at  the  interface.  Also  for  the  remaining  elements  one 
obtains  the  coefficients  in  terms  of  the  previously  calculated  values  of 

bK*  - 

For  the  Kth  element 
p(K)(x)  =  i»K  +  bK  X  +  CK 


Requiring  the  derivatives  to  be  equal  at  the  interface: 
P(K)(0)  =  bK  =  2cj<_1  L  +  bK_]  =  p(K'V)(L) 

Requiring  the  functions  to  be  continuous: 

P(K)(0)  =  aK  =  pK 
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P^K  ]\l)  =  pk_*,  4  bK_1  L  +  cK-1  ,_2  =  p 


Solving  for  aK,  bR  and  c.K  one  finds; 


'  PK 


bK  =  “bK-l  +  2 


(pK  "  Pk-1^ 


_  bK-l  ,  W+l  "  3pK  +  2pK-1^ 

Ck  _____  +  - p - “ 


A  similar  analysis  for  the  yelocity  gives 
v(l)(x)  =  ^  +  51  X2 


where 


51  =  V1 


n,  = 


~v3  +  4v2 


-  3Vn 


2L 


_  V3  -  2v2  +  V1 


2L2 


and 


v(K)(x)  =  +  \  x  +  5k  x2 


where 


?K  =  VK 
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nK  =  '  Vi  +  2 


<VK  -  Vl> 


V  ,  (vfc 


„  _  K-1  L  VVK+1 
+ - 


-  3v»  +  2vl  ' 


'K 

IT 


'K-1  ■ 


Also,  the  temperature  distribution  is  handled  in  the  same  way: 

T^(x)  =  6-j  +  0-jX  + 

where 

61 

*1 

Y1 
and 

\x)  =  6k  +  V  +  V2 
where 

6K 

yK 


=  Ti, 


=  -  *K-1  +  2 


5K-1 


<tk  -  W 


(Ti 


K+l 


L 

"k 


IT  +  °T 
*  *  “  K  V 


IT 


=  Tn 


-T3  +  4T2  -  31-j 
2L 


T3  -  2T2  +  T1 
2L2 


One  now  applies  the  state,  nomentor,  continuity,  and  energy  o^ua- 
tions  to  ojcain  relationships  fcr  t,.e  pressure,  teno  ratrrv 

vjlocity,  a..c!  density  at  eac.i  nod  or  interface  prirvl. 


frs 


-  Total  Mass  for  Kth  element 

P^K^tot  -  Total  Momentum  for  Kth  element 
(Kl 

IP  Of  ”  ^"ota^  Potential  energy  for  Kth  element 
K^tqt  “  ^°tal  Kinetic  energy  for  Kth  element 


From  previous  work,  a^,  b^,  and  are  known  as  functions  of 
p.j,  pg...pK  except  for  K  =  1,  which  is  a  function  of  p^  also. 

One  obtains  n  differential  equations,  i.e.,  an  equation  for  each 
element.  Similarily,  one  obtains  n  equation's  from  the  momentum  and 
energy  relationships. 

Momentum  equation  is: 

d 

a  r  tot 

dt  =  "  A  ^PK+1  VK+1  “PKVK2  +  PK+1  "  PK^ 
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where 


P^TOT  =  A  Ax)  pK(x)  dx 


PK  is  eliminated  by  using  the  state  equation: 


PK  -  RpK  Tk 


After  expanding  and  integrating 


•  /  i  ,  _  lA 

aK  (ckL  +  \  2  +  h  3  ] 


+  6k  (ck  L2  +  \  3  +  CK  +, 

•  /  L3  I  4  l5\ 

CK  T+  nK  T+  CK 


^aK  L  +  bK  2  +CK  A 


•  ,  L2  L3  l4, 

nK  ^aK  T  +  bK  X  +  CK  f 

•  ,a  L3  .  L4  ,  „  LA 

5K  ^  K  3  bK  4  CK  5  ^  “ 


^PK+1  VK+1  PKVK  +  ^pk+1  TX+1  “  ^p((  V 


The  energy  equation  is: 


d  (u(k)tot  +  K(K)- 

dt 


"  A  °  ^p;c  vk  "  pic+i  vk+!i  tk+i^ 


4 


)  +  0<K>  -  *<K> 


+  1  (P|<  VK3  ‘  PK+1  vl+] 

where 

p<K)(x)  <v(K)(x))2  d* 

0 

fL 

m(K)  _  »  /  D\  0(K)  t(K) 

l  T0T  -  A  (a  R)  p  (x)  T  (x)  ax 

o 

-  rate  heat  enters  an  element 

-  work  being  done  within  an  element 

aK  L  +  2?KnKT+  ^K5K  +  nVl  +  ?nK5KT 

+  ?2kT)/2  +  {a-R)  <L6k+  hT+  ykT)]  + 

^;2KT+  2?KnKT  +  +  nVV 

+  Z\hT  +  ?2kT^/2  +  (a“R)  ^6KT+  6KT+  yKT'^ 

+  CK  ^c2KT+  2cKr'KT+  ^KCK  +  n2K^5/5 

+  2nKCKT  +  c2KT'/2  +  ^6kT+  6kV+  ykT^ 

+  [cKaKL  +  UKt>K  1-  nKaK)4r+  UKcK  +  nKbK  +  eKa«) 

+  (nk/C|/  +  +  s./Cix-r-] 


^2  f  j  ^ 

*  \  ^?KaKT+  ^?KbK  +  ,,KaK^T+  (?KCK  +  nKbK 

+  EKaK>T  +  <\eK  +  EKbK>T  +  VKT] 

+  EK  ^K8kV  +  ^KbK  +  nKaK^T+  WCK  +  nKbK 

+  EKaK>T  +  ("KCK  +  5KbK)T'  +  «KcKT’  f 
6|/  C(ot-R)  (a^L  +  b^-  +  cKy  )] 

+  6k  [(a-R)  (dK^-+  bK^-+  CK^)] 

+  YK  C(«-R)  («Ky  +  bK^+  CKy)] 

=  «  (pkvkTk  -  P(<+1  vK+]  Tk+1)  +  £•  (pKvK3  -  pK+1  v3+1) 

< 

(0(K)  -  fl(K))/A 

With  the  derived  equations,  the  computer  can  be  used  to  generate  the 
coefficients  of  the  derivatives  and  then  to  solve  for  the  derivatives. 


